/******************************************************************************/
* Authors: 	A.M. Birkenbach, M.-Y. Lee, & M.D. Smith
* Date: 	08/02/2023
* Project: 	Counterfactual Modeling of Multispecies Fisheries Outcomes Under 
*			Market-Based Regulation
* Purpose: 	Produce graphs/other tables in paper
* Inputs:	actuals_for_comparison.dta 
*			monthly_summary_validation1.csv
*			monthly_summary_validation2.csv
*			monthly_summary_counterfactual_closemultB.csv
*			econ_202*.csv
/******************************************************************************/

clear all
gl input ""
gl working ""
gl graphs ""
gl output ""


********************************************************************************
*define program "dataprep"
cap program drop dataprep
program define dataprep
	*create monthly count of nofish decisions by replicate and model
	preserve
	keep if spstock2=="nofish"
	/* 	r_ = landings*prices
		p_ = landings*lag prices */
	keep fy_* model replicate gearcat trips l_ c_ r_ spstock2
	duplicates drop
	rename trips t_
	rename spstock2 sp
	qui tempfile nofish
	qui save `nofish', replace
	restore
	append using `nofish'
	collapse (sum) l_ c_ r_ t_ /*q_*/, by(sp *month *year gearcat replicate model)
	*generate distribution across replicates
	collapse (p50) l_50=l_ c_50=c_ r_50=r_ t_50=t_ ///
		(p25) l_25=l_ c_25=c_ r_25=r_ t_25=t_ ///
		(p75) l_75=l_ c_75=c_ r_75=r_ t_75=t_, ///
		by(sp *month *year gearcat model)
	*merge in actuals
	merge m:1 fy_year fy_month gearcat sp using "$input\actuals_for_comparison.dta", keep(1 3)
	assert _merge==3 if sp~="nofish"
	replace month=fy_month+4 if fy_month<=8 & month==.
	replace month=fy_month-8 if fy_month>8 & month==.
	replace year=fy_year if month>=5 & year==.
	replace year=fy_year-1 if month<=4 & year==.
	drop _merge
	gen monthyear=mdy(month,1,year)
	format monthyear %tdMonCCYY
	*label vars
	rename choice t_actual
	rename qtykept l_actual
	rename revenue r_actual
	foreach var of varlist *_actual {
		label var `var' "Actual (Catch Shares)"
	}
	label var monthyear "Month"
	foreach i in "t" "l" "c" "r" {
		forvalues j=25(25)75 {
			label var `i'_`j' "Predicted (p`j')"
		}
	}
	renpfix l_ landings_
	renpfix c_ catch_
	renpfix r_ revenues_ // quota/DAS charges NOT subtracted out
	renpfix t_ days_fished_
	order model gearcat sp sppname_new monthyear
	gen monthyearm=mofd(monthyear)
	format monthyearm %tm
	label var monthyearm "Month"
	replace sppname_new="No Fish" if sp=="nofish"
end

********************************************************************************
*compile simulation results
cd "$input"
foreach name in "validation1" "validation2" "counterfactual_closemultB" {
	import delim "$input\monthly_summary_`name'.csv", varn(1) clear
	tempfile t`name'
	save `t`name'', replace
}
clear
gen source=""
foreach name in "validation1" "validation2" "counterfactual_closemultB" {
	append using `t`name''
	replace source="`name'" if source==""
}
tab model
drop v1 *rev_total* *choice_prev_fish
foreach var of varlist month year {
    rename `var' fy_`var'
}
replace model=subinstr(model,"counterfactual","counterfactual_closemult",1) if strpos(source,"mult")~=0
replace model="validation" if inlist(model,"validation1","validation2") 
replace model=subinstr(model,"mult","multB",1) if strpos(source,"multB")~=0
tab model
drop source
reshape long l_ c_ r_ q_, i(*month *year spstock2 gearcat replicate model trips)j(sp)string
*note: reshape creates gearcat_sp combos that shouldn't exist. delete.
replace spstock2=lower(spstock2)
gen t_ = trips if spstock2==sp
bysort gearcat sp: egen totaltrips=total(t_)
sum totaltrips
assert `r(min)'==0
tab sp gearcat if totaltrips==0
drop if totaltrips==0
compress
keep if	inlist(model,"counterfactual_closemultB_pre_coefsnc2","validation_pre_coefsnc2"), clear
dataprep
replace model=substr(model,1,1)
compress
reshape wide *_25 *_50 *_75, i(gearcat sp monthyear* month year *_actual) j(model)string
egen group=group(gearcat sp*)
tsset group monthyearm, monthly
foreach var of varlist *_actual {
	label var `var' "Actual (CS)"
}
foreach i in "landings" "revenues" "days_fished" {
	label var `i'_50c "Predicted (DAS)" // c=counterfactual
	label var `i'_50v "Predicted (CS)"  // v=validation
	label var `i'_25c "25th percentile (DAS)"
	label var `i'_25v "25th percentile (CS)"
	label var `i'_75c "75th percentile"
	label var `i'_75v "75th percentile"
}
gen landings_CS_DAS = landings_50v - landings_50c
gen revenues_CS_DAS = revenues_50v - revenues_50c
bysort gearcat: egen totalrevdiff=total(revenues_CS_DAS)
gen gf=inlist(sp,"americanlobster","monkfish","other","redsilveroffshorehake","seascallop")==0 & ///
	inlist(sp,"skates","spinydogfish","squidmackerelbutterfishherring","summerflounder","nofish")==0
bysort gearcat gf: egen totalrevdiffgf=total(revenues_CS_DAS)
format totalrevdiff* %12.0g
save "$working\tempgraphdata.dta", replace



foreach name in "validation" "validation2" "closemultB" {
	cd "$input"
	local allfiles : dir . files "econ_202*.csv"
	clear
	foreach f in `allfiles' {
		import delim using "`f'", varn(1) clear
		loc g=strtoname("`f'")
		gen trips=1
		collapse (sum) c_* l_* r_* trips, by(hullnum spstock2 gearcat m) 
		tempfile `g'
		save `g', replace
	}
	foreach f in `allfiles' {
		loc g=strtoname("`f'")
		append using `g'.dta
		erase `g'.dta
	}
	save "$input\vesselrev_`name'.dta", replace
}
use "$input\vesselrev_validation.dta", clear
append using "$input\vesselrev_validation2.dta"
append using "$input\vesselrev_closemultB.dta"
replace model="counterfactual_closemult" if model==""
replace m=subinstr(m,"counterfactual","counterfactual_closemult",1) if strpos(model,"mult")~=0
tab m
drop model
rename m model
*save dataset of trip counts by hullnum
preserve
collapse (sum) trips, by(hullnum gearcat model spstock2) 
replace trips=trips/100
save "$input\
.dta", replace
restore
collapse (sum) r_* l_* c_*, by(hullnum gearcat model) // already summed over replicates, now summing over time
*to take the average over all 100 replicates, divide by 100
foreach var of varlist c_* r_* l_* {
	replace `var'=`var'/100
}
reshape long r_ l_ c_, i(model gearcat hullnum)j(spstock)string
*note: reshape creates gearcat_sp combos that should not exist. delete.
bysort gearcat spstock: egen sumrev=total(r_)
drop if sumrev==0
tab spstock gearcat
drop sumrev
save "$input\vesselrev.dta", replace
erase "$input\vesselrev_closemultB.dta"
erase "$input\vesselrev_validation.dta"
erase "$input\vesselrev_validation2.dta"


********************************************************************************
*Paired t-tests
cap log close
use "$input\vesselrev.dta", clear
encode spstock, gen(sp)
assert r_~=.
qui levelsof model, loc(models)
qui levelsof hullnum, loc(hulls)
qui log using "$input\Gini_byvessel.txt", replace text // excluded from DATAVERSE for confidentiality
foreach m in `models' {
	foreach h in `hulls' {
		preserve
		qui keep if model=="`m'" & hullnum=="`h'"
		qui levelsof gearcat, loc(g) clean
		qui replace r_=0.0000001 if r_==0
		qui ginidesc r_, by(sp) mat(gini)
		di "`m', `g', `h',", gini[1,1]	
		restore
	}
}
qui cap log close
insheet using "$input\Gini_byvessel.txt", clear nonames comma
egen rowmiss=rowmiss(v1-v4)
keep if rowmiss==0
drop rowmiss
rename v1 model
rename v2 gearcat
rename v3 hullnum
rename v4 gini
replace gini=. if gini==0
keep if inlist(model,"validation_pre_coefsnc2","counterfactual_closemult_pre_coefsnc2")
replace model=subinstr(model,"_"," ",.)
replace model=word(model,1)
compress
reshape wide gini, i(gearcat hullnum)j(model)string
gen ginidiff=ginivalidation-ginicounterfactual
cap erase "$input\pairedttests.txt"

ttest ginivalidation==ginicounterfactual
matrix ttest= (r(N_1), r(mu_1), r(mu_2), r(mu_1)-r(mu_2), r(p))
matrix rownames ttest="ALL"
matrix colnames ttest= pairs meancs meandas mean_difference pscore
mat2txt, matrix(ttest) sav("$input\pairedttests.txt") append

ttest ginivalidation==ginicounterfactual if gearcat=="GILLNETS"
matrix ttest= (r(N_1), r(mu_1), r(mu_2), r(mu_1)-r(mu_2), r(t), r(p))
matrix rownames ttest="GILLNETS"
matrix colnames ttest= pairs meancs meandas mean_difference tstat pscore
mat2txt, matrix(ttest) sav("$input\pairedttests.txt") append

ttest ginivalidation==ginicounterfactual if gearcat=="TRAWL"
matrix ttest= (r(N_1), r(mu_1), r(mu_2), r(mu_1)-r(mu_2), r(t), r(p))
matrix rownames ttest="TRAWL"
matrix colnames ttest= pairs meancs meandas mean_difference tstat pscore
mat2txt, matrix(ttest) sav("$input\pairedttests.txt") append

********************************************************************************
*Gini of total revenues over vessels in each gearcat (i.e., concentration of revenues across vessels)
cap log close
use "$input\vesselrev.dta", clear
collapse (sum) r_, by(model gearcat hullnum model)
encode hullnum, gen(hull)
qui levelsof model, loc(models)
qui levelsof gearcat, loc(gears)
qui log using "$input\Gini_bygearcat.txt", replace text
foreach g in `gears' {
	foreach m in `models' {
		preserve
		qui keep if model=="`m'" & gearcat=="`g'"
		qui replace r_=0.0000001 if r_==0
		qui ginidesc r_, by(hull) mat(gini)
		di "`g', `m',", gini[1,1]	
		restore
	}
}
qui cap log close
import delim using "$input\Gini_bygearcat.txt", delim(",") clear
drop if inlist(v1,"TRAWL","GILLNETS")==0
compress
rename v1 gearcat
rename v2 model
rename v3 gini
keep if strpos(model,"nc2")~=0
drop if strpos(model,"closeown")~=0



cap log close
qui log using "$output\4__GRAPHS_log.txt", replace text

********************************************************************************
*VALIDATION REGRESSIONS
use "$input\simresults.dta" if strpos(model,"validation")~=0, clear
dataprep
encode sp, gen(sp2)
replace model=subinstr(model,"validation_pre_coefs","Model ",1)
levelsof model, loc(models)
levelsof gearcat, loc(gears)
foreach m in `models' {
    foreach g in `gears' {
		qui reg landings_actual landings_50 i.sp2 if model=="`m'" & gearcat=="`g'", ro
		loc r2new=round(`e(r2)',.01)
		di "`m', `g', `r2new'"
	}
}

********************************************************************************
*Figures 1 & 2: Validation simulations of monthly landings
*Figures 3 & 4: Counterfactual simulations of monthly landings
*Figures 8 & 9: Fishing days by target stock under CS vs. DAS
*validation and counterfactual line graphs
use "$working\tempgraphdata.dta", clear
levelsof gearcat, l(gears)
foreach gear in `gears' {
	preserve
	keep if gearcat=="`gear'"
	drop if sppname=="No Fish"
	sort sp monthyear
	*Landings: Validation
	graph twoway (rarea landings_75v landings_25v monthyear, lc(gs11) fc(gs11)) ///
		(line landings_actual monthyear, lc(black)) ///
		(line landings_50v monthyear, lc(gs8) lw(medium)), ///
		by(sppname_new, style(compact) yr note("") graphr(c(white))) subtitle(, size(small)) ///
		scheme(s1manual) ytitle("Landings (lbs)") xlabel(,angle(45)) tlabel(#7) tmtick(##4) ///
		ylabel(,angle(horizontal) labs(small) glc(gs14)) 
	qui gr export "$graphs\\`gear'landingsvalidationALL.png", replace width(1600) height(900)
	*Landings: Counterfactual
	graph twoway (rarea landings_75c landings_25c monthyear, lc(gs11) fc(gs11)) ///
		(line landings_50c monthyear, lc(gs8)) ///
		(rarea landings_75v landings_25v monthyear, lc(gs4) fc(gs4)) ///
		(line landings_50v monthyear, lc(black) lw(medium)), ///
		by(sppname_new, style(compact) yr note("") graphr(c(white))) subtitle(, size(small)) ///
		scheme(s1manual) ytitle("Landings (lbs)") xlabel(,angle(45)) tlabel(#7) tmtick(##4) ///
		ylabel(,angle(horizontal) labs(small) glc(gs14)) 
	qui gr export "$graphs\\`gear'landingscounterfactualALL.png", replace width(1600) height(900)
	*Fishing days: Counterfactual
	drop if sppname=="No Fish" // exit not consistently defined as no-fish across CS and DAS, so leave no-fish days out
	sort sppname_new monthyear
	graph twoway line days_fished_50c days_fished_50v monthyear, ///
		by(sppname_new, yr note("") style(compact) graphr(c(white)) subtitle(, size(small))) ///
		lc(black gs10) lw(medium medium) scheme(s1manual) ytitle("Fishing Vessel-Days per Month") ///
		xlabel(,angle(45) grid) tlabel(#7) tmtick(##4) ylabel(,angle(horizontal))
	qui gr export "$graphs\\`gear'daysfished.png", replace width(1600) height(900)			
	restore
}

********************************************************************************
*Figure 5: Cumulative distribution of landings over the year under CS vs. DAS
use "$working\tempgraphdata.dta", clear
collapse (sum) landings_50*, by(sp gf month sppname_new)
gen fy_month=month+8 if month<5
replace fy_month=month-4 if month>=5
drop if sp=="nofish"
foreach var of varlist landings_50* {
	bysort sp (fy_month): gen cumsum`var'=sum(`var')
	bysort sp (fy_month): egen total`var'=total(`var')
	gen pctotal`var'=cumsum`var'/total`var'
	drop total`var'
}
gen uniform=1/12*fy_month
label var pctotallandings_50v "Predicted (CS)"
label var pctotallandings_50c "Predicted (DAS)"
label var fy_month "Month of Fishing Year"
label var uniform "Uniform"
gen gf_inv=1-gf // to get gf to show up first
graph twoway line pctotallandings_50c pctotallandings_50v* uniform fy_month /*if gf==1*/, ///
	by(gf_inv sppname_new, yr note("") holes(15) graphr(c(white))) subtitle(, size(small)) scheme(s1manual) ///
	ytitle("Avg. Portion of Total Yearly Landings") ///
	ylabel(,angle(horizontal)) cmissing(n n n) ///
	xlabel(1(2)12) lc(black gs11 gs9) lw(medthick medthick vthin) ///
	lp(solid solid shortdash) legend(row(1) region(lstyle(none)))
qui gr export "$graphs\cumulqty.png", ///
	replace width(1400) height(1300)
	
********************************************************************************
*Figure 6. : Simulated groundfish revenue differences (CS-DAS)
use "$working\tempgraphdata.dta", clear
collapse (sum) revenues_CS_DAS, by(monthyear sppname_new sp gf)
graph twoway area revenues_CS_DAS monthyear if gf==1 & inlist(sp,"yellowtailfloundersnema","winterfloundergom")==0, ///
	by(sppname_new, style(compact) note("") graphr(c(white))) subtitle(, size(medsmall)) ///
	scheme(s1manual) xlabel(,angle(45)) tlabel(#7) tmtick(##4) ///
	ylabel(,angle(horizontal) labs(small) glc(gs14)) 
qui gr export "$graphs\revdiffsgf.png", replace width(1600) height(1000)

********************************************************************************
*Figure 7: Decomposition of simulated revenue differences (CS-DAS) into volume and price components
use "$working\tempgraphdata.dta", clear
drop if sp=="nofish"
keep gearcat sppname_new fy* revenues_50* landings_50* gf
renpfix revenues_50 revenues_
renpfix landings_50 landings_
foreach model in "v" "c" {
	gen price_`model'=revenues_`model'/landings_`model'
}
gen price_c_landings_v=price_c*landings_v
gen price_v_landings_c=price_v*landings_c
collapse (sum) revenues* price*landings*, by(gearcat gf /*sppname_new*/)
gen vi=.5*(revenues_v - price_v_landings_c + price_c_landings_v - revenues_c)
gen pi=.5*(revenues_v - price_c_landings_v + price_v_landings_c - revenues_c)
ds, has(type numeric)
format price* revenues* %12.0g
assert abs((vi+pi)-(revenues_v-revenues_c))<5
foreach var of varlist vi pi {
	replace `var'=`var'/1000000
}
order gearcat gf vi pi
cap label drop gf
label define gf 0 "Non-Groundfish" 1 "Groundfish"
label values gf gf
graph bar vi pi, over(gf) over(gearcat) ///
	legend(label(1 "Volume component") label(2 "Price component") region(c(none))) ///
	scheme(s1manual) ylabel(,angle(horizontal)) plotregion(margin(zero)) ///
	bar(1, color(gs11)) bar(2, color(gs5)) ytitle("Contribution to revenue differences (CS - DAS) ($ millions)", size(small))
qui gr export "$graphs\bennet.png", replace width(1400) height(1000)

********************************************************************************
*Figures 10 & 11: Revenue composition under CS vs. DAS
*Groundfish
use "$input\vesselrev.dta", clear
foreach st in "ccgom" "gb" "gom" "snema" {
	replace spstock=subinstr(spstock,"`st'","",1)
}
tab spstock
collapse (sum) r_, by(model gearcat spstock)
replace model="CS" if model=="validation_pre_coefsnc2"
replace model="DAS" if model=="counterfactual_closemult_pre_coefsnc2"
keep if inlist(model,"CS","DAS")
gen gf=inlist(sp,"americanlobster","monkfish","other","redsilveroffshorehake","seascallop")==0 & ///
	inlist(spstock,"skates","spinydogfish","squidmackerelbutterfishherring","summerflounder")==0
replace spstock=proper(spstock)
replace spstock=subinstr(spstock,"flounder"," Fl.",1)
tab spstock
replace spstock="White Hake" if spstock=="Whitehake"
replace spstock="SMB-Herring" if spstock=="Squidmackerelbutterfishherring"
replace spstock="Spiny Dogfish" if spstock=="Spinydogfish"
replace spstock="Sea Scallop" if spstock=="Seascallop"
replace spstock="RSO Hake" if spstock=="Redsilveroffshorehake"
replace spstock="American Lobster" if spstock=="Americanlobster"
replace spstock="Am. Plaice" if spstock=="Americanplaice Fl."

cap erase "bar1g.gph"
cap erase "bar1t.gph"
graph hbar (sum) r_ if gearcat=="GILLNETS" & gf==1, over(spstock) over(model) ///
	asy per stack blabel(bar, position(center) format(%9.0f) color(white)) ///
	graphr(c(white)) plotregion(lc(none)) scheme(s1manual) ytitle("") ylab(,labs(0) tl(0)) ///
	legend(size(small) span rows(1) region(lc(white))) pcycle(10) ///
	title("GILLNETS") saving(bar1g)
graph hbar (sum) r_ if gearcat=="TRAWL" & gf==1, over(spstock) over(model) ///
	asy per stack blabel(bar, position(center) format(%9.0f) color(white)) ///
	graphr(c(white)) plotregion(lc(none)) scheme(s1manual) ytitle("Revenue composition (%)") ///
	legend(size(small) span rows(1) region(lc(white))) pcycle(10) ///
	title("TRAWL") saving(bar1t)
gr combine bar1g.gph bar1t.gph, xcom ycom col(1) ///
	scheme(s1manual) ysize(10) xsize(16)
qui gr export "$graphs\revhbar1.png", as(png) width(2000) replace

*Groundfish and non-groundfish
replace spstock="Groundfish" if gf==1
encode spstock, gen(spstock2)
label values spstock2 .
replace spstock2=1 if spstock=="Groundfish"
replace spstock2=2 if spstock=="American Lobster"
labmask spstock2, values(spstock)
qui cap erase "bar2g.gph"
qui cap erase "bar2t.gph"
graph hbar (sum) r_ if gearcat=="GILLNETS", over(spstock2) over(model) ///
	asy per stack blabel(bar, position(center) format(%9.0f) color(white)) ///
	graphr(c(white)) plotregion(lc(none)) scheme(s1manual) ytitle("") ylab(,labs(0) tl(0)) ///
	legend(size(small) span rows(1) region(lc(white))) pcycle(10) ///
	title("GILLNETS") saving(bar2g)
graph hbar (sum) r_ if gearcat=="TRAWL", over(spstock2) over(model) ///
	asy per stack blabel(bar, position(center) format(%9.0f) color(white)) ///
	graphr(c(white)) plotregion(lc(none)) scheme(s1manual) ytitle("Revenue composition (%)") ///
	legend(size(small) span rows(1) region(lc(white))) pcycle(10) ///
	title("TRAWL") saving(bar2t)
gr combine bar2g.gph bar2t.gph, xcom ycom col(1) ///
	scheme(s1manual) ysize(10) xsize(15)
qui gr export "$graphs\revhbar2.png", as(png) width(2000) replace

********************************************************************************

qui log close


*end of do-file